#本模块兼具求解pde笔记的作用
#FiPy（FVM），GetFem++（FEM），FEniCS/DOLFIN（FEM）
#Theoretical and Numerical Background--from fipy官网-Theoretical and Numerical Background
#https://www.ctcms.nist.gov/fipy/documentation/numerical/index.html
#   有限元法使用FV单元的特征函数（即等于1的函数）作为加权函数，可以得到与FVM完全等价的方程组[17]。
#   类似地，用有限体积法对方程进行离散化可以简化为笛卡尔网格上的有限体积法。

#一些常见的参数代表含义
'''
rank(参数)--变量（CellVariable，FaceVariable等）的维度
            --（the rank (number of dimensions) of each element of this Variable. Default: 0
CellVariable(对象)
    0.参数：mesh name  value rank elementshape unit 
    1.value属性：numpy.ndarray对象
    2.mesh属性：mesh=CellVariable挂载的mesh对象
    3.()操作：即CellVariable对象的__call__(points=None, order=0, nearestCellIDs=None)属性
                插值返回 points点的数值
                point: tuple或者list或者ndarray 表示求值的点坐标 ( (0., 1.1, 1.2), (0., 1., 1.) ) 表示三个点(0,1);(1.1,1);(1.2,1)
                order:插值阶数，默认为0（返回最近点）. order=1:好像是线性插值
    4.print操作：直接返回 CellVariable.value的结果
    5.copy()操作：
    6.faceGrad属性：Return \nabla \phi as a rank-1 FaceVariable using differencing for the normal direction(second-order gradient).
    
    7.[]操作：返回一个能标记位置的东西--<class 'fipy.variables.unaryOperatorVariable._UnaryOperatorVariable.<locals>.unOp'>
            如CellVariable[0]，其mesh会取一个index的限制--mesh=UniformGrid2D(dx=1.0, nx=3, dy=1.0, ny=2))[index] == 0.0
                print(CellVariable[0])会得到 x坐标
                > 操作：如 CellVariable[0] > 1 返回同一种对象，此时有 [index] >= 1

    7.constrain属性：phi.constrain(valueLeft, mesh.facesLeft)
                mesh.facesLeft返回的对象作用和[]操作差不多，都能标记位置

FaceVariable(对象)
    以mesh为 UniformGrid2D(dx=1.0, nx=3, dy=1.0, ny=2)为例
    0. 大部分与CellVariable相同
        nx=3,  ny=2 对应17条边界(所有边界)
    1. divergence属性：与CellVariable的faceGrad相对应

mesh(对象)--网格
    以 UniformGrid2D(dx=1.0, nx=3, dy=1.0, ny=2)为例
    1. cellCenters属性：CellVariable(对象),记录每个网格中心的坐标值.
        [[0.5 1.5 2.5 0.5 1.5 2.5]
         [0.5 0.5 0.5 1.5 1.5 1.5]] shape为(2,6) 从左下到右上，先x值，再y值
    2 facesLeft属性：类似FaceVariable的对象
            其mesh会取一个index ：mesh=UniformGrid2D(dx=1.0, nx=3, dy=1.0, ny=2))[index] == 0.0
 
variable(对象):
    变量
    1. value属性：ndarray 值
    2. ()操作:  返回value属性
ExponentialConvectionTerm

'''

#先考虑Fipy求解
from fipy import Variable, FaceVariable, CellVariable, Grid1D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer
from fipy.tools import numerix

nx = 50
dx = 1.
mesh = Grid1D(nx=nx, dx=dx)
D = 1.
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
valueLeft = 1
valueRight = 0

phi.constrain(valueRight, mesh.facesRight)
phi.constrain(valueLeft, mesh.facesLeft)
eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D)
timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100

phiAnalytical = CellVariable(name="analytical value", mesh=mesh)
if __name__ == '__main__':
    viewer = Viewer(vars=(phi, phiAnalytical), datamin=0., datamax=1.)
    viewer.plot()
x = mesh.cellCenters[0]
t = timeStepDuration * steps
try:
    from scipy.special import erf
    phiAnalytical.setValue(1 - erf(x / (2 * numerix.sqrt(D * t))))
except ImportError:
    print("The SciPy library is not available to test the solution to \
the transient diffusion equation")

for step in range(steps):
    eqX.solve(var=phi, dt=timeStepDuration)
    if __name__ == '__main__':
        viewer.plot()
